本来想做绕轴旋转,绕轴旋转是在垂直于轴向量的空间里旋转,但是n维空间里与某个向量垂直的空间为n-1维,而旋转只在二维空间里有定义。所以这里就改成了,向任意方向旋转。
计算单位向量 X = ( x 1 , x 2 , ⋯ , x n ) X=(x_1,x_2,\cdots,x_n) X=(x1,x2,⋯,xn)旋转到单位向量 V = ( v 1 , v 2 , ⋯ , v n ) V=(v_1,v_2,\cdots,v_n) V=(v1,v2,⋯,vn)的旋转矩阵 R R R。
旋转坐标系
任意的旋转都可以看作绕着一个轴,在某个平面上的旋转。
不失一般性,假定向量
V
=
(
v
1
,
v
2
,
⋯
,
v
n
)
V=(v_1,v_2,\cdots,v_n)
V=(v1,v2,⋯,vn)在
v
1
×
v
2
v_1\times v_2
v1×v2张成的平面上旋转,用矩阵乘法可表示为:
V
′
=
R
1
V
V'=R_1V
V′=R1V其中旋转矩阵
R
1
R_1
R1定义为
R
1
=
(
cos
α
−
sin
α
sin
α
cos
α
1
⋱
1
)
R_1=\left( \begin{array}{ccc} \cos\alpha& -\sin\alpha& & & \\ \sin\alpha& \cos\alpha& & & \\ & &1 & & \\ & & & \ddots & \\ & & & & 1 \\ \end{array} \right)
R1=⎝⎜⎜⎜⎜⎛cosαsinα−sinαcosα1⋱1⎠⎟⎟⎟⎟⎞
其中
cos
α
=
v
2
v
1
2
+
v
2
2
\cos\alpha=\frac{v_2}{\sqrt{v_1^2+v_2^2}}
cosα=v12+v22v2,
sin
α
=
v
1
v
1
2
+
v
2
2
\sin\alpha=\frac{v_1}{\sqrt{v_1^2+v_2^2}}
sinα=v12+v22v1。
由此我们将向量
V
V
V旋转到与
(
v
1
,
0
,
⋯
,
0
)
(v_1,0,\cdots,0)
(v1,0,⋯,0)正交。得
V
′
=
(
0
,
v
1
2
+
v
2
2
,
v
3
,
v
4
,
⋯
,
v
n
)
V'=(0,\sqrt{v_1^2+v_2^2}, v_3,v_4,\cdots,v_n)
V′=(0,v12+v22,v3,v4,⋯,vn)。
以此类推,将
V
V
V旋转到
(
0
,
⋯
,
0
,
1
)
(0,\cdots,0,1)
(0,⋯,0,1)需要
n
−
1
n-1
n−1次旋转。将这些旋转变换组合为一个旋转矩阵
R
0
=
R
n
−
1
R
n
−
2
⋯
R
1
R_0=R_{n-1}R_{n-2}\cdots R_1
R0=Rn−1Rn−2⋯R1(注意,乘法不满足交换律,顺序不能颠倒)。顺便我们可以计算出这个旋转的逆矩阵
R
0
−
1
=
R
1
⊤
⋯
R
n
−
2
⊤
R
n
−
1
⊤
R_0^{-1}=R_1^\top\cdots R_{n-2}^\top R_{n-1}^\top
R0−1=R1⊤⋯Rn−2⊤Rn−1⊤。
旋转向量
现在我们要做的是将 R 0 X R_0X R0X旋转到坐标轴 ( 0 , ⋯ , 0 , 1 ) (0,\cdots,0,1) (0,⋯,0,1),并求旋转矩阵。做法与上面相同。最后可以求得,在旋转后的坐标系下,旋转矩阵为 R X R_X RX: R 0 V = R X R 0 X R_0V=R_XR_0X R0V=RXR0X
虽然
R
R
R和
R
X
R_X
RX表示的旋转角度相同,但是这两个旋转所在的坐标系是不同的,还需要将
R
X
R_X
RX的坐标轴旋转回去:
R
=
R
0
−
1
R
X
R
0
R=R_0^{-1}R_XR_0
R=R0−1RXR0
代码实现
import numpy as np
def Rot_map(V):
assert len(V.shape) == 1
assert np.linalg.norm(V) - 1 < 1e-8
n_dim = V.shape[0]
Rot = np.eye(n_dim)
Rot_inv = np.eye(n_dim)
for rotate in range(n_dim-1):
rot_mat = np.eye(n_dim)
rot_norm = np.sqrt(V[rotate]**2 + V[rotate+1]**2)
cos_theta = V[rotate+1]/rot_norm
sin_theta = V[rotate]/rot_norm
rot_mat[rotate,rotate] = cos_theta
rot_mat[rotate,rotate+1] = - sin_theta
rot_mat[rotate+1,rotate] = sin_theta
rot_mat[rotate+1,rotate+1] = cos_theta
V = np.dot(rot_mat, V)
Rot = np.dot(rot_mat, Rot)
Rot_inv = np.dot(Rot_inv,rot_mat.transpose())
return Rot, Rot_inv
n_dim = 512
X = np.random.rand(n_dim)
X = X / np.linalg.norm(X)
V = np.random.rand(n_dim)
V = V / np.linalg.norm(V)
R_0, R_0_inv = Rot_map(V)
R_X, _ = Rot_map(np.dot(R_0, X))
R = np.dot(np.dot(R_0_inv, R_X), R_0)
assert np.linalg.norm(np.dot(R,X) - V) < 1e-8
扩展
这个方法可以用来在正交于某个向量的空间中随机采样。例如,在 n n n维空间中,采样与 V = ( v 1 , v 2 , ⋯ , v n ) V=(v_1,v_2,\cdots,v_n) V=(v1,v2,⋯,vn)正交的向量,则代码如下:
import numpy as np
def Rot_map(V):
assert len(V.shape) == 1
assert np.linalg.norm(V) - 1 < 1e-8
z_dim = V.shape[0]
Rot = np.eye(z_dim)
Rot_inv = np.eye(z_dim)
for rotate in range(z_dim-1):
rot_mat = np.eye(z_dim)
rot_norm = np.sqrt(V[rotate]**2 + V[rotate+1]**2)
cos_theta = V[rotate+1]/rot_norm
sin_theta = V[rotate]/rot_norm
rot_mat[rotate,rotate] = cos_theta
rot_mat[rotate,rotate+1] = - sin_theta
rot_mat[rotate+1,rotate] = sin_theta
rot_mat[rotate+1,rotate+1] = cos_theta
V = np.dot(rot_mat, V)
Rot = np.dot(rot_mat, Rot)
Rot_inv = np.dot(Rot_inv,rot_mat.transpose())
return Rot, Rot_inv
z_dim = 512
n_samples = 10000
V = np.random.rand(z_dim)
V = V / np.linalg.norm(V)
R_0, R_0_inv = Rot_map(V)
samples = np.random.rand(z_dim-1, n_samples)
samples = np.concatenate((samples, np.zeros((1, n_samples))))
samples = np.dot(R_0_inv, samples)
assert np.mean(np.abs(np.dot(V, samples))) < 1e-8
进一步,我们采样与一组向量 V s Vs Vs正交的向量:
import numpy as np
def gs(X, row_vecs=True, norm = True):
'''
https://gist.github.com/iizukak/1287876/edad3c337844fac34f7e56ec09f9cb27d4907cc7
'''
if not row_vecs:
X = X.T
Y = X[0:1,:].copy()
for i in range(1, X.shape[0]):
proj = np.diag((X[i,:].dot(Y.T)/np.linalg.norm(Y,axis=1)**2).flat).dot(Y)
Y = np.vstack((Y, X[i,:] - proj.sum(0)))
if norm:
Y = np.diag(1/np.linalg.norm(Y,axis=1)).dot(Y)
if row_vecs:
return Y
else:
return Y.T
def Rot_map(V):
assert len(V.shape) == 1
assert np.linalg.norm(V) - 1 < 1e-8
z_dim = V.shape[0]
Rot = np.eye(z_dim)
Rot_inv = np.eye(z_dim)
for rotate in range(z_dim-1):
rot_mat = np.eye(z_dim)
rot_norm = np.sqrt(V[rotate]**2 + V[rotate+1]**2)
cos_theta = V[rotate+1]/rot_norm
sin_theta = V[rotate]/rot_norm
rot_mat[rotate,rotate] = cos_theta
rot_mat[rotate,rotate+1] = - sin_theta
rot_mat[rotate+1,rotate] = sin_theta
rot_mat[rotate+1,rotate+1] = cos_theta
V = np.dot(rot_mat, V)
Rot = np.dot(rot_mat, Rot)
Rot_inv = np.dot(Rot_inv,rot_mat.transpose())
return Rot, Rot_inv
z_dim = 10 # number of dimensions
n_vectors= 2 # number of Vs
assert n_vectors < z_dim
n_samples = 10 # number of samples to orthogonal Vs
Vs = np.random.rand(n_vectors, z_dim)
# Gram-Schmidt Orthogonization
orth_Vs = gs(Vs)
orth_Vs = orth_Vs.T
Rots = np.eye(z_dim)
Rot_invs = np.eye(z_dim)
for idx in range(n_vectors):
V = orth_Vs[:, idx]
R_0 = np.eye(z_dim)
R_0_inv = np.eye(z_dim)
R_0[0:z_dim-idx,0:z_dim-idx], R_0_inv[0:z_dim-idx,0:z_dim-idx] = Rot_map(V[range(z_dim-idx)])
orth_Vs = R_0.dot(orth_Vs)
Rots = np.dot(R_0, Rots)
Rot_invs = np.dot(Rot_invs, R_0_inv)
samples = np.random.rand(z_dim-n_vectors, n_samples)
samples = np.concatenate((samples, np.zeros((n_vectors, n_samples))))
samples = np.dot(Rot_invs, samples)
assert np.mean(np.abs(np.dot(Vs, samples))) < 1e-8
用矩阵旋转太慢了,上述采样操作可以用矩阵分解解决(在任意向量的正交空间中采样):
import numpy as np
z_dim = 10 # number of dimensions
n_vectors= 2 # number of Vs
assert n_vectors < z_dim
n_samples = 20 # number of samples to orthogonal Vs
Vs = np.random.rand(n_vectors, z_dim)
Vs = np.concatenate((Vs,np.zeros((z_dim-n_vectors, z_dim))))
U, s, V = np.linalg.svd(Vs)
directions = np.eye(z_dim)
directions[range(n_vectors),range(n_vectors)] = 0
directions = np.dot(np.dot(U, directions), V)
samples = np.random.rand(z_dim, n_samples)
samples = np.dot(directions.T, samples)
assert np.mean(np.abs(np.dot(Vs, samples))) < 1e-8